Both the time series show a linear upward trend. Visually, the Paperback seems to have some regular patters (seasonality) approximately equal to 3 days.
autoplot(books, facets = T)
There doesn’t seem to be any linear correlation between the two series.
lattice::xyplot(Paperback~Hardcover, as.data.frame(books), type=c('p','smooth'))
The PACF plots do show that for Paperback, there is a 3-order autocorrelation in the signal which is significant. For Hardcover, there is a 1st and 2nd order significance.
ggPacf(books[,'Paperback'])
ggPacf(books[,'Hardcover'])
Here are three settings - a=0.9, a=0.5, and a=0.01. As alpha increases, so does the uncertainty of the prediction since ses will look back further in time. Though at alpha = 0.01, the point estimate seems quite low.
ses(books[,'Paperback'], initial = 'simple', alpha = .9) %>% forecast() %>% autoplot()
ses(books[,'Paperback'], initial = 'simple', alpha = .5) %>% forecast() %>% autoplot()
ses(books[,'Paperback'], initial = 'simple', alpha = .01) %>% forecast() %>% autoplot()
We can see a typical curve as seen during parameter tuning. The SSE is minimum at an alpha value of about 0.2.
sse_list <- c()
a_list <- seq(0.001, 0.999, length.out = 20)
for (a_sel in a_list) {
ses(books[,'Paperback'], initial = 'simple', alpha = a_sel)$model$SSE -> sse
sse_list <- c(sse_list, sse)
}
plot(a_list, sse_list, type='b')
points(a_list[which.min(sse_list)], sse_list[which.min(sse_list)], col='red', pch=20)
Alpha selected is 0.215. The point estimate for this forecast seem better than the #2 results. Prediction interval widths are about similar. RMSE is smaller for the auto-alpha selection.
ses(books[,'Paperback'], initial = 'simple', alpha = NULL)$model$par['alpha']
ses(books[,'Paperback'], initial = 'simple', alpha = NULL)%>% forecast(h=4)%>% autoplot()
bind_rows(ses(books[,'Paperback'], initial = 'simple', alpha = NULL) %>%
forecast(h=4) %>% accuracy() %>% sweep::sw_tidy(),
ses(books[,'Paperback'], initial = 'simple', alpha = 0.01) %>%
forecast(h=4) %>% accuracy() %>% sweep::sw_tidy()) %>%
knitr::kable(format = 'latex',caption = '1st line: Auto-alpha. 2nd line: a=0.01')
Setting the initial to optimal reduces RMSE by ~1 unit, and MASE by ~0.02 units.
ses(books[,'Paperback'], initial = 'optimal', alpha = NULL) -> fit_optimal
ses(books[,'Paperback'], initial = 'simple', alpha = NULL) -> fit_simple
bind_rows(fit_optimal %>% forecast(h=4) %>% accuracy() %>% sweep::sw_tidy(),
fit_simple %>% forecast(h=4) %>% accuracy() %>% sweep::sw_tidy()) %>%
knitr::kable(format = 'latex',
caption = '1st line: Optimal initial. 2nd line: Simple initial')
This series does better with a higher value of alpha.
ses(books[,'Hardcover'], initial = 'simple', alpha = .9) %>% forecast() %>% autoplot()
ses(books[,'Hardcover'], initial = 'simple', alpha = .5) %>% forecast() %>% autoplot()
ses(books[,'Hardcover'], initial = 'simple', alpha = .01) %>% forecast() %>% autoplot()
We can see a typical curve as seen during parameter tuning. The SSE is minimum at an alpha value of about 0.35.
sse_list <- c()
a_list <- seq(0.001, 0.999, length.out = 20)
for (a_sel in a_list) {
ses(books[,'Hardcover'], initial = 'simple', alpha = a_sel)$model$SSE -> sse
sse_list <- c(sse_list, sse)
}
plot(a_list, sse_list, type='b')
points(a_list[which.min(sse_list)], sse_list[which.min(sse_list)], col='red', pch=20)
Alpha selected is 0.347. The point estimate for this forecast seem better than previous results. Prediction interval widths are about similar. RMSE is smaller for the auto-alpha selection.
ses(books[,'Hardcover'], initial = 'simple', alpha = NULL)$model$par['alpha']
ses(books[,'Hardcover'], initial = 'simple', alpha = NULL)%>% forecast(h=4)%>% autoplot()
bind_rows(ses(books[,'Hardcover'], initial = 'simple', alpha = NULL) %>%
forecast(h=4) %>% accuracy() %>% sweep::sw_tidy(),
ses(books[,'Hardcover'], initial = 'simple', alpha = 0.9) %>%
forecast(h=4) %>% accuracy() %>% sweep::sw_tidy()) %>%
knitr::kable(format = 'latex',caption = '1st line: Auto-alpha. 2nd line: a=0.01')
Setting the initial to optimal reduces RMSE by ~2 units, but MASE increased by ~0.01 units.
ses(books[,'Hardcover'], initial = 'optimal', alpha = NULL) -> fit_optimal
ses(books[,'Hardcover'], initial = 'simple', alpha = NULL) -> fit_simple
bind_rows(fit_optimal %>% forecast(h=4) %>% accuracy() %>% sweep::sw_tidy(),
fit_simple %>% forecast(h=4) %>% accuracy() %>% sweep::sw_tidy()) %>%
knitr::kable(format = 'latex',
caption = '1st line: Optimal initial. 2nd line: Simple initial')
holt(books[,'Paperback']) %>% forecast(h=4) %>% autoplot()
holt(books[,'Hardcover']) %>% forecast(h=4) %>% autoplot()
The holt method definitely has a lower SSE than the sse method since holt can account for the upward trend in both the timeseries. ses unfortunately, cannot account for this trend.
ses(books[,'Paperback']) %>% residuals() %>% .^2 %>% sum()
ses(books[,'Hardcover']) %>% residuals() %>% .^2 %>% sum()
holt(books[,'Paperback']) %>% residuals() %>% .^2 %>% sum()
holt(books[,'Hardcover']) %>% residuals() %>% .^2 %>% sum()
Again, the forecasts created by holt have the right upward trend as would be expected from the
Using the forecast function:
holt(books[,'Paperback']) %>% forecast(h=1)
holt(books[,'Hardcover']) %>% forecast(h=1)
autoplot(ukcars)
stl_fit <- stl(ukcars, s.window = 'periodic')
autoplot(stl_fit)
sadjusted <- seasadj(stl_fit)
holtfit <- holt(sadjusted, h = 8, damped = T, exponential = F)
seasonaladjustments <- ukcars-sadjusted
holtfit_w_seasonal <- holtfit$mean + seasonaladjustments[2:9]
holtfit_w_seasonal
Parameters of the fit are here. RMSE is 25.15986.
summary(holtfit)
holtfit_additive <- holt(sadjusted, h = 8, damped = F, exponential = F)
seasonaladjustments <- ukcars-sadjusted
holtfit_additive_w_seasonal <- holtfit_additive$mean + seasonaladjustments[2:9]
holtfit_additive_w_seasonal
Parameters of the fit are here. RMSE is 25.26.
summary(holtfit_additive)
ETS selects an A-Ad-A model - Additive error, Additive damped seasonal and Additive trend component.
ets(ukcars, model = 'ZZZ',damped = T) %>% summary()
The ETS gives marginally worse results.
Given how close the RMSE values are we expect very similar forecasts. And we can see this in the plots. They are virtually indistinguishable.
holtfit$mean <- holtfit$mean + seasonaladjustments[2:9]
holtfit$upper <- holtfit$upper + seasonaladjustments[2:9]
holtfit$lower <- holtfit$lower + seasonaladjustments[2:9]
holtfit$x <- holtfit$x + seasonaladjustments
ets(ukcars, model = 'ZZZ',damped = T) %>% forecast() %>% autoplot()
holtfit %>% autoplot()
autoplot(visitors)
hw(visitors, seasonal = 'm', damped = F) %>%
forecast(h=12*2) %>% autoplot()
This is because the seasonality keeps increasing over time. We can visually see this if we keep the seasonality constant (using decompose). Look at the residuals - the magnitude keeps increasing over time.
Multiplicative seasonality allows for it to increase over time.
decompose(visitors) %>% plot
As expected, the exponential method overpredicts the point estimates than the damped forecast. But, the exponential method seems to have a smaller prediction interval.
hw(visitors, seasonal = 'm', damped = T) %>%
forecast(h=12*2) %>% autoplot()
hw(visitors, seasonal = 'm', exponential = T) %>%
forecast(h=12*2) %>% autoplot()
The RMSEs are fairly similar for all three models. Based on RMSE I would pick the damped Holt Winters model. MASE for this model is the lowest too.
bind_rows(
hw(visitors, seasonal = 'm') %>% forecast(h=12*2) %>%
accuracy() %>% sweep::sw_tidy(),
hw(visitors, seasonal = 'm', damped = T) %>% forecast(h=12*2) %>%
accuracy() %>% sweep::sw_tidy(),
hw(visitors, seasonal = 'm', exponential = T) %>% forecast(h=12*2) %>%
accuracy() %>% sweep::sw_tidy()
) %>% mutate(.rownames = c('hw_m','hw_m_d','hw_m_e')) %>%
rename(model = .rownames)
hw_m <- hw(visitors, seasonal = 'm') %>% forecast(h=12*2)
hw_m %>% autoplot()
ets_fit <- ets(visitors) %>% forecast(h=12*2)
ets_fit %>% autoplot()
ets_box_fit <- ets(visitors, lambda = 'auto') %>% forecast(h=12*2)
ets_box_fit %>% autoplot()
snaive_box_fit <- snaive(visitors, lambda = 'auto') %>% forecast(h=12*2)
snaive_box_fit %>% autoplot()
This code performs the needed actions. The stl function removes the trend from the signal very well. The residuals look fairly random visually.
After adjusting for seasonality, we can see that an ETS model is an A-Ad-N model: no seasonality with a linear damped trend. Now, although we have an additive trend component, the beta coeef is 1e-4, so practically, the forecast is flat as we can see in the plot.
BoxCox.lambda(x = visitors)
visitors_boxed <- BoxCox(x = visitors, lambda = 0.2775249)
visitors_stl <- stl(visitors_boxed, s.window = 'periodic')
plot(visitors_stl)
visitors_sadj <- seasadj(visitors_stl)
ets_sadj <- ets(visitors_sadj, model = 'ZZZ', damped = T)
ets_sadj %>% summary()
ets_sadj %>% forecast() %>% autoplot()
The residuals tell us this:
checkresiduals(hw_m)
checkresiduals(ets_fit)
checkresiduals(ets_box_fit)
checkresiduals(snaive_box_fit)
checkresiduals(ets_sadj)
Putting everything together, we can see that:
Given this information and the residuals, I would probably pick the ETS model. It has almost no autocorrelation and seems to fit the data the best.
hw_m = hw_m %>% forecast(h = 12*2) %>% sweep::sw_tidy() %>%
pull('Point Forecast') %>% ts(frequency = 12, start=c(2005,5))
ets_fit = ets_fit %>% forecast(h = 12*2) %>% sweep::sw_tidy() %>%
pull('Point Forecast') %>% ts(frequency = 12, start=c(2005,5))
ets_box_fit = ets_box_fit %>% forecast(h = 12*2) %>% sweep::sw_tidy() %>%
pull('Point Forecast') %>% ts(frequency = 12, start=c(2005,5))
snaive_box_fit = snaive_box_fit %>% forecast(h = 12*2) %>% sweep::sw_tidy() %>%
pull('Point Forecast') %>% ts(frequency = 12, start=c(2005,5))
ets_reseasoned = ets_sadj %>% forecast(h = 12*2) %>% sweep::sw_tidy() %>%
pull('Point Forecast') %>% ts(frequency = 12, start=c(2005,5))
ets_reseasoned = (ets_reseasoned + visitors_stl$time.series %>%
tail(24) %>% .[,'seasonal'] %>%
as.numeric()) %>% InvBoxCox(lambda = 0.2775249)
ts.union(visitors, hw_m, ets_fit, ets_box_fit, snaive_box_fit, ets_reseasoned) %>%
autoplot(size=2) + xlim(c(2000,2008))
y <- ts(data = numeric(100))
e <- rnorm(n = 100, mean = 0, sd = 1)
for(i in 2:100)
y[i] <- 0.6*y[i-1] + e[i]
plot(y)
phi1 <- seq(-1, 1, by = 0.5)
phi1
[1] -1.0 -0.5 0.0 0.5 1.0
y <- ts(data = numeric(100))
e <- rnorm(n = 100, mean = 0, sd = 1)
result <- matrix(nrow = 100,ncol = length(phi1))
for(p in seq_along(phi1)){
for(i in 2:100)
y[i] <- phi1[p]*y[i-1] + e[i]
result[,p] <- y
}
plot (result[,1], type='l', col='gray')
lines(result[,2], type='l', col='red')
lines(result[,3], type='l', col='blue')
lines(result[,4], type='l', col='green')
lines(result[,5], type='l', col='orange')
Over the range of values of phi_1 from -1 to 1:
y <- ts(data = numeric(100))
e <- rnorm(n = 100, mean = 0, sd = 1)
for (i in 2:100){
y[i] <- 0.6*e[i-1] + e[i]
}
plot(y)
There doesn’t seem to any concrete relationship between theta and the result.
y <- ts(numeric(100))
e <- rnorm(100)
for(i in 2:100)
y[i] <- 0.6 * y[i - 1] + 0.6 * e[i - 1] + e[i]
plot(y)
y <- ts(numeric(100))
e <- rnorm(100)
for(i in 3:100)
y[i] <- 0.8 * y[i - 1] + 0.3 * y[i - 2] + e[i]
plot(y)
Graphed above. While the ARMA(1,1) gives a stationary series, the AR(2) is clearly non-stationary.
kpss.test(wmurders)
p-value smaller than printed p-value
KPSS Test for Level Stationarity
data: wmurders
KPSS Level = 1.2017, Truncation lag parameter = 1, p-value = 0.01
The KPSS test shows us that we must reject the H0: data are stationary. Differencing the time series by 1:
ggtsdisplay(diff(wmurders,lag = 1))
kpss.test(diff(wmurders,1))
KPSS Test for Level Stationarity
data: diff(wmurders, 1)
KPSS Level = 0.58729, Truncation lag parameter = 1, p-value = 0.02379
Trying one more differencing scheme. Double differencing - Differencing the difference by 1:
ggtsdisplay(diff(diff(wmurders,lag = 1)))
Thus, I can guess that an ARIMA(0, 2, 2) might fit the data well.
Yes, I think we should include a constant. With a constant value included, with a d=1 or d=2, a long term forecast can follow a line or quadratic curve. The time series seems to be downward trending, so a constant will help. If c = 0, with d = 0 or d = 1, we won’t get a downward trend. (With c = 0 and d = 2, we can still get a linear trend).
(1 - B)^2 y_t = c + (1 + theta1 x B + theta2 x B^2) x e_t
fit1 <- Arima(wmurders, order = c(0, 2, 2))
summary(fit1)
Series: wmurders
ARIMA(0,2,2)
Coefficients:
ma1 ma2
-1.0181 0.1470
s.e. 0.1220 0.1156
sigma^2 estimated as 0.04702: log likelihood=6.03
AIC=-6.06 AICc=-5.57 BIC=-0.15
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set -0.0113461 0.2088162 0.1525773 -0.2403396 4.331729 0.9382785 -0.05094066
checkresiduals(fit1)
Ljung-Box test
data: Residuals from ARIMA(0,2,2)
Q* = 11.764, df = 8, p-value = 0.1621
Model df: 2. Total lags used: 10
The Ljung-Box test shows a p-val > 0.05, so there is little chance of autocorrelation in the residuals. Visually looking at the ACF plot, we can see some there is periodic elements in the residuals. The model seems adequate.
fit1 %>% forecast(h = 3)
Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
2005 2.480525 2.202620 2.758430 2.055506 2.905544
2006 2.374890 1.985422 2.764359 1.779250 2.970531
2007 2.269256 1.772305 2.766206 1.509235 3.029276
fit1 %>% forecast(h = 3) %>% autoplot()
If we allow stepwise:
auto.arima(wmurders, seasonal = F, allowdrift = T)
Series: wmurders
ARIMA(1,2,1)
Coefficients:
ar1 ma1
-0.2434 -0.8261
s.e. 0.1553 0.1143
sigma^2 estimated as 0.04632: log likelihood=6.44
AIC=-6.88 AICc=-6.39 BIC=-0.97
If we force a thorough investigation:
auto.arima(wmurders, seasonal = T, stepwise = F, allowdrift = T)
Series: wmurders
ARIMA(0,2,3)
Coefficients:
ma1 ma2 ma3
-1.0154 0.4324 -0.3217
s.e. 0.1282 0.2278 0.1737
sigma^2 estimated as 0.04475: log likelihood=7.77
AIC=-7.54 AICc=-6.7 BIC=0.35
The ARIMA(0,2,3) has the lowest AICc value (-6.7) against my selection of ARIMA(0,2,2) with an AICc of -5.7.
I’m going to log transform the data
Very strong seasonality - 4, 8, 12, 16 seen in the components. Reducing strength of the contributions as expected.
Strong seasonal components at lag=4. Perhaps a strong non-seasonal component at lag=5.
After adjusting for seasonality, and differencing the data once to make it stationary, the ACF and PACF plots tell me a significant lag of 1 exists. It’s hard for me to pick of p=1 or q=1 based on this plot alone.
Perhaps I might choose a ARIMA(1,0,0)(0,1,1)[4] or a ARIMA(0,0,1)(0,1,1)[4]. But I can’t pick without more analysis.
auto.arima(laustourists, stepwise = F)
Series: laustourists
ARIMA(1,0,0)(0,1,1)[4] with drift
Coefficients:
ar1 sma1 drift
0.4154 -0.9043 0.0128
s.e. 0.1387 0.2711 0.0011
sigma^2 estimated as 0.004541: log likelihood=54.52
AIC=-101.05 AICc=-100.02 BIC=-93.91
It’s similar to what I picked, and I tend to agree with it’s selection.
(1 - phi1 x B)(1 - PHI1 x B^4)(1 - B^4)y_t = (1 + THETA1 x B^4)e_t
Yes, there is an increasing seasonality. Looks like an inverse sqrt will work.
BoxCox.lambda(usmelec)
[1] -0.4772402
The data are not stationary. A difference of 1 is needed to make it stationary according to the KPSS test.
diff(lusmelec,1) %>% kpss.test()
p-value greater than printed p-value
KPSS Test for Level Stationarity
data: .
KPSS Level = 0.022613, Truncation lag parameter = 4, p-value = 0.1
Visually, the signal looks stationary.
Zooming in a bit more, I think I can see seasonality which a diff=1 hasn’t got rid of.
If we look at the ACF, PACF plots, we can see two patterns: * a seasonal pattern at 12, 24, 36 * another seasonal pattern at 3, 6, 9, …
Looks like we have a complex dual-seasonal pattern even after differencing.
The PACF plot shows lags at 12, 24, 36 which are seasonal lags. Perhaps an AR(3) term for seasonal is needed. For non-seasonal component, considering the drop in ACF values, the most significant PACF value is for lag = 1. Perhaps an ARIMA(1,0,0)(0,1,3)[12] is a good guess.
It looks like ARIMA(1,0,1)(1,1,1)[12] has the lowest AIC of -4240.98, after some investigation.
Arima(lusmelec, order = c(1,0,0), seasonal = list(order=c(0,1,3), period=12))
Series: lusmelec
ARIMA(1,0,0)(0,1,3)[12]
Coefficients:
ar1 sma1 sma2 sma3
0.9843 -0.8726 -0.0957 0.1128
s.e. 0.0106 0.0512 0.0581 0.0502
sigma^2 estimated as 4.358e-06: log likelihood=2098.67
AIC=-4187.35 AICc=-4187.21 BIC=-4166.89
Arima(lusmelec, order = c(1,0,0), seasonal = list(order=c(0,1,2), period=12))
Series: lusmelec
ARIMA(1,0,0)(0,1,2)[12]
Coefficients:
ar1 sma1 sma2
0.9858 -0.8526 -0.0223
s.e. 0.0100 0.0551 0.0531
sigma^2 estimated as 4.398e-06: log likelihood=2096.22
AIC=-4184.43 AICc=-4184.34 BIC=-4168.07
Arima(lusmelec, order = c(1,0,0), seasonal = list(order=c(1,1,2), period=12))
Series: lusmelec
ARIMA(1,0,0)(1,1,2)[12]
Coefficients:
ar1 sar1 sma1 sma2
0.9864 -0.4831 -0.3578 -0.4591
s.e. 0.0099 0.3726 0.3533 0.3042
sigma^2 estimated as 4.401e-06: log likelihood=2096.53
AIC=-4183.06 AICc=-4182.92 BIC=-4162.6
Arima(lusmelec, order = c(1,0,1), seasonal = list(order=c(0,1,3), period=12))
Series: lusmelec
ARIMA(1,0,1)(0,1,3)[12]
Coefficients:
ar1 ma1 sma1 sma2 sma3
0.9965 -0.5114 -0.8071 -0.1279 0.0895
s.e. 0.0037 0.0645 0.0512 0.0577 0.0508
sigma^2 estimated as 3.881e-06: log likelihood=2126.49
AIC=-4240.98 AICc=-4240.79 BIC=-4216.44
Arima(lusmelec, order = c(1,0,1), seasonal = list(order=c(1,1,1), period=12))
Series: lusmelec
ARIMA(1,0,1)(1,1,1)[12]
Coefficients:
ar1 ma1 sar1 sma1
0.9971 -0.5191 0.0702 -0.8720
s.e. 0.0033 0.0645 0.0568 0.0309
sigma^2 estimated as 3.903e-06: log likelihood=2124.75
AIC=-4239.51 AICc=-4239.37 BIC=-4219.05
fit_selected <- Arima(lusmelec, order = c(1,0,1), seasonal = list(order=c(1,1,3), period=12))
fit_selected %>% summary()
Series: lusmelec
ARIMA(1,0,1)(1,1,3)[12]
Coefficients:
ar1 ma1 sar1 sma1 sma2 sma3
0.9959 -0.5069 0.3740 -1.1789 0.1630 0.1226
s.e. 0.0042 0.0651 0.3026 0.3007 0.2565 0.0506
sigma^2 estimated as 3.881e-06: log likelihood=2127.04
AIC=-4240.07 AICc=-4239.81 BIC=-4211.43
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.0001088471 0.001930585 0.001511984 0.00564225 0.07789407 0.5754233 0.1284949
fit_selected %>% checkresiduals()
Ljung-Box test
data: Residuals from ARIMA(1,0,1)(1,1,3)[12]
Q* = 54.46, df = 18, p-value = 1.556e-05
Model df: 6. Total lags used: 24
Residuals show a left skew and the ACF plot shows significant lags at 1,2,14 etc. Trying out another model: an ARIMA(2,0,2)(1,1,3)[12] improved the AIC to -4264.
fit_selected <- Arima(lusmelec, order = c(2,0,2), seasonal = list(order=c(1,1,3), period=12))
fit_selected %>% summary()
Series: lusmelec
ARIMA(2,0,2)(1,1,3)[12]
Coefficients:
ar1 ar2 ma1 ma2 sar1 sma1 sma2 sma3
1.3027 -0.3035 -0.7186 -0.0789 0.4392 -1.2694 0.2476 0.1115
s.e. 0.1826 0.1822 0.1872 0.1259 0.3036 0.3025 0.2690 0.0489
sigma^2 estimated as 3.656e-06: log likelihood=2141.09
AIC=-4264.18 AICc=-4263.77 BIC=-4227.36
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set 3.731734e-05 0.001869441 0.001448626 0.001974489 0.07462978 0.5513111 0.02330224
fit_selected %>% checkresiduals()
Ljung-Box test
data: Residuals from ARIMA(2,0,2)(1,1,3)[12]
Q* = 29.996, df = 16, p-value = 0.01802
Model df: 8. Total lags used: 24
What does an auto.arima say? ARIMA(1,1,1)(2,1,1)[12]… but the AICc is worse than the model I chose.
fit_selected <- auto.arima(lusmelec, stepwise = F, approximation = F)
fit_selected <- auto.arima(lusmelec, stepwise = F, approximation = F)
fit_selected %>% summary()
Series: lusmelec
ARIMA(1,1,1)(2,1,1)[12]
Coefficients:
ar1 ma1 sar1 sar2 sma1
0.4002 -0.8295 0.0269 -0.1016 -0.8485
s.e. 0.0652 0.0385 0.0579 0.0553 0.0366
sigma^2 estimated as 3.653e-06: log likelihood=2135.33
AIC=-4258.65 AICc=-4258.46 BIC=-4234.12
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set -9.026359e-05 0.001873122 0.001423716 -0.004660559 0.07332098 0.5418307 0.007791069
fit_selected %>% checkresiduals()
Ljung-Box test
data: Residuals from ARIMA(1,1,1)(2,1,1)[12]
Q* = 27.628, df = 19, p-value = 0.09086
Model df: 5. Total lags used: 24
The forecasted model works quite well!! The green line shows the real data for the past few years, while the red line is my forecasts. I am quite pleased.
Probably a year or two years out. Because forecasting assumes the underlying data generating process remains unchanged. If this changes, then the model is useless.